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Abstract 

The osmotic equation of state for the athermal bond fluctuation model on the simple cubic lattice is 
obtained from extensive Monte Carlo simulations. For short macromolecules (chain length N=20) we 
study the influence of various choices for the chain stiffness on the equation of state. Three techniques 
are applied and compared in order to critically assess their efficiency and accuracy: the "repulsive 
wall" method, the thermodynamic integration method (which rests on the feasibility of simulations 
in the grand canonical ensemble), and the recently advocated sedimentation equilibrium method, 
which records the density profile in an external (e.g. gravitation-like) field and infers, via a local 
density approximation, the equation of state from the hydrostatic equilibrium condition. We confirm 
the conclusion that the latter technique is far more efficient than the repulsive wall method, but we 
find that the thermodynamic integration method is similarly efficient as the sedimentation equilibrium 
method. For very stiff chains the onset of nematic order enforces the formation of an isotropic-nematic 
interface in the sedimentation equilibrium method leading to strong rounding effects and deviations 
from the true equation of state in the transition regime. 
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I. INTRODUCTION 



Understanding the equation of state of macromolecules in solution has been a longstanding 
problem of polymer science, which is both important as a fundamental problem in the statis- 
tical mechanics of soft matter-'2i>^'^i^i^>^ and relevant for various applications of polymers. The 
interplay of excluded volume effects^'^'^'^, solvent quality and variable chain stiffness-'-^ already 
is very difficult to describe for macromolecules in very dilute solution^, withstanding an analytic 
solution and making the application of computer simulation methods^-iifi necessary. The exper- 
imental situation on measurements of the equation of state of polymer solutions is addressed 
in chapter 5 of Ref. and chapter 3 of Ref. 11. Considering semidilute and concentrated so- 
lutions, the enthalpic and/or entropic interactions among the polymer chains create nontrivial 
correlations in the structure of the solutions involving many chains, and phase transitions such 
as phase separation under bad solvent conditions into a dilute polymer solution and a concen- 
trated one may occur-i2i»i^'i^»i^. Alternatively, under good solvent conditions one may observe 
for semifiexible or for stiff chains a transition from an isotropic to a nematic solution when 
the polymer concentration increases^'i^'i^'i^ii^. Of particular interest is the situation when ne- 
matic ordering and the tendency to phase separation under bad solvent conditions compete^^. 
Preliminary computer simulation studies of a corresponding model2ii2^ gave only rather rough 
information on the phase behavior, and it was concluded that the accuracy with which the 
osmotic equation of state could be determined needs to be improved. 

Thus, there are many reasons why establishing methods that allow the reliable estimation of 
the equation of state for various coarse-grained lattice models for polymers is of interest. While 
in an off-lattice model the estimation of the pressure tensor from the virial theorem in principle 
is straightforward^^, it often is not possible to study very large systems with the desirable 
accuracy^, and hence simulations of lattice models still have their place^ii^. However, for lattice 
models different techniques to calculate the osmotic pressure of polymers in solution must be 
sought in order to derive their equation of state. The standard method, particularly valuable 
for dilute solutions and/or not too large chain lengths, calculates the chemical potential from 
the insertion probability of a test chain2ii2^ and the osmotic pressure then follows from standard 
thermodynamic integration^i^iS^i^. As is well known, insertion of a polymer chain into a volume 
containing already other chains is very difficult due to very small acceptance probabilities, and 
hence requires advanced Monte Carlo methods2^i22i>2E to be practically feasible. 

As an alternative method Dickman^^'^^ proposed the repulsive wall thermodynamic integra- 
tion (RWTI) method, which remains applicable also for very dense systems since it works in the 
canonic {MVT) ensemble, where the number of chains M in the box remains fixed. However, 
this method requires substantial simulation effort, since for each state point in the bulk several 
simulation runs with increasing wall-monomer repulsion and subsequent thermodynamic inte- 
gration need to be performed. Moreover the RWTI method may suffer from important finite 
size effects^. 

More recently, it was proposed^^ to extract the osmotic equation of state from Monte Carlo 
simulations where the equilibrium monomer and center-of-mass concentration profiles of lattice 
polymers in a gravitation-like potential are computed. From these concentration profiles, the 
equation of state can be inferred if a local density approximation like in hydrostatic equilibrium 
is invoked. This method has been broadly applied to study the sedimentation equilibrium of 
colloidal dispersions containing spherical^, rod-like^^ or disk-like^ particles. More recently, 
successful applications to colloid-polymer mixturesSZ and solutions of block copolymers'^ or 
binary polymer solutions'^^ have been made as well. 

Despite these successes, this "sedimentation equilibrium" (SE) method also may have draw- 
backs, when other large length scales appear in the system, that compete with the characteris- 
tic "sedimentation length". This length that scales inversely with the "gravitational constant" 
characterizing the gravitation-like potential and hence can be made arbitrarily large^^, but the 
linear dimension of the simulated system in the 2;-direction in which this gravitation-like force 



acts must then be huge, also. A well-known case, where even the real gravitation potential on 
earth substantially disturbs the equation of state is a fluid very close to the gas-liquid critical 
point^, since there the correlation length of density fluctuations diverges, and critical fluctua- 
tions undisturbed by gravity can occur in the x, y-directions perpendicular to the gravitational 
force only™. Other cases, apart from critical points of second-order transitions, where large 
length scales arise that are potentially disturbed by the gravitation-like potential are associated 
with the formation of thick wetting layers at walls, for instance. 

While Addison at al.— did test the SE method against the RWTI method for a few cases 
varying the solvent quality from good solvents to theta solvents, flnding good agreement, they 
considered only fully flexible chains. Thus, we take up this problem but rather consider semi- 
flexible chains as well: the possible occurrence of nematically ordered wetting layers at the wall 
where the density is largest^l or even the occurrence of the isotropic to nematic transition in 
the bulk solution of the semiflexible or stiff polymers may complicate matters. The aim of 
the present work is to carefully test the SE method against both the RWTI method and the 
standard thermodynamic integration method in the grand canonical fiVT ensemble [TlfiVT 
method, denoting the chemical potential), and to assess by detailed comparisons both the 
efficiency and the accuracy of these methods. As has been explained above, there is need for 
accurate simulation methods for the computation of the equation of state of polymer solutions 
in the context of various interesting problems. 

The plan for the remainder of this paper is as follows. In Sec. [IT] we briefly review the 
theoretical basis for the RWTI, Tl^VT and SE methods, while Sec. [111] describes the simulated 
model and mentions the types of Monte Carlo moves used. Sec. [IV] presents comparisons of 
the three methods for fully flexible chains over a wide range of densities, while Sec. IVl presents 
our results for variable chain stiffness. Sec. IVII contains our conclusions and gives an outlook 
to future work. 



II. METHODS TO CALCULATE THE OSMOTIC PRESSURE FOR LATTICE MOD- 
ELS OF POLYMER SOLUTIONS 

A. Thermodynamic integration in the grand canonical ensemble (TlfiVT method) 

We consider a system of chains of length in a simulation box of volume V with periodic 
boundary conditions (typically the box has a cubic shape, V = L^, L being the linear dimension) 
at given temperature T and chemical potential /i of the chains. The density of polymer chains 
in the system, p = Af/V, with Af the average number of chains contained in the simulation box 
then is an output of the simulation. 

Utilizing the standard relations in the canonical ensemble, fi = {d{F/V)/dp)Ty and p = 
— {dF / dV)T,Mi where F is the free energy of the system and p the pressure, one easily derives 
that 

p 

P = pp — J l^{p')dp' + const, (1) 



where the integration constant in Eq. ([1]) can be flxed by reference to the low density limit, 
where the system behaves like an ideal gas of chains, 

pV = AfkeT, vr = p/keT = p. (2) 

Denoting then the chemical potential of the ideal gas of chains as pid, and deflning p'^^ = 
p — pid the excess chemical potential per chain, Eqs. ([I]),© imply^ 

p 

7r = p(l + /i-)- I p'^p'W (3) 





In practice the integral in Eq. ([3]) is discretized, so that the reduced osmotic pressure tTj at 
chain density pi is obtained from the recursion relation 

TT, ^ 7r,_i + (1 + /ir)p. - (1 + ^A-l)p^~l - (/^r + ^^T-l){p^ - A-i)/2. (4) 

The disadvantages of the method are clearly obvious from Eq. (jlj): a large number of state points 
{/ij, T, V} needs to be studied with small differences between /ij and and hence pi and pj_i, 
so that the discretization error going from Eq. ^ to Eq. (jlj) is negligible; and an efficient grand 
canonical simulation method is needed, so pi is sampled with sufficient accuracy. For not too 
long chains and not too high densities, however, the method is indeed practically useful, and 
since periodic boundary conditions are used, the system for large V always is homogeneous, 
and the analysis is not hampered by any interfacial effects, which are present inevitably in the 
other techniques described below due to their explicit use of walls. 



B. The repulsive wall thermodynamic integration (RWTI) method 

This method can be implemented both in the grand canonical pVT ensemble and in the 
canonical AfVT ensemble. While the latter choice has the advantage that no chain insertions 
are necessary and hence the method also works for long chains and at high densities, it suffers 
from rather large finite size effects^^. We have used here both the canonical and the grand 
canonical version of the method. 

One considers now a box of linear dimensions L x L x H , with periodic boundary conditions 
in X and y directions only, while a hard wall is placed both at z = and z = H. Moreover, one 
introduces a repulsive potential of strength Swaii which acts in the first layer adjacent to z = 
and in the last layer available to the monomers, z = H. As discussed in2Si2ii22^ the osmotic 
pressure is then obtained from the fraction of sites, 02(A), occupied by monomers in the layers 
adjacent to both walls, ^ = and z = H, 

1 

TT = / y I 1 , A = exp{-e^,all/kBT). (5) 



Again the integration in Eq. ([5]) is discretized, and 20 different values of A turned out to be 
sufficient to obtain reliable results. 



C. The sedimentation equilibrium (SE) method 

This method utilizes the canonical ensemble N'VT, and one also considers a box of linear 
dimensions LxLxH with periodic boundary conditions in x and y directions only, while again 
hard walls are used in z-direction at z = and sX z = H. An external potential is applied not 
at the walls but everywhere in the system^ 

Ue^tern^\{z) = -UigZ = ~'^9^- (6) 

Here m is the mass of a monomer, g is the acceleration due to the gravity-like potential, 
a is the lattice spacing, and A^ is a dimensionless constant that characterizes the strength 
of this gravitational potential, Xg = amg/ksT. It is also useful to introduce characteristic 
gravitational lengths ^m, icm 

^m = 0'/Xg, ^cm = ^m/N. (7) 

As shown by Addison et al.^^, for large z the density profile of an ideal gas of monomers at the 
lattice would follow the standard barometric formula, Pm{z) oc exp{—mgz/ ksT) = exp{—z/^rn), 



while for an ideal gas of polymer chains the density profile of monomer units Pm{z) = Np{z) oc 

While the variation of the density profile for large z, where the system is very dilute and 
the ideal gas behavior holds, hence is trivially known, and this knowledge is an important 
consistency check of the method, for smaller z the density profile is nontrivial. But from this 
profile the osmotic equation of state can be estimated, when one invokes the local density 
approximation, such that the equation of hydrostatic equilibrium holds— 

^ = -Nmgpiz), (8) 
p{z) being the local osmotic pressure at altitude z. Integration of Eq.® yields 



IT 



oo 

{z)=p{z)/kBT = i;^ j p{z')dz'. (9) 



Thus one can record both p{z) (and the associated monomer profile Pm{z)) and 7r(z) simul- 
taneously, and eliminating z from these relations one obtains the desired equation of state 7r(p) 
(or n{pm), respectively). 

Noting that in the local density approximation we expect that the density distribution of 
the monomer units satisfies Pm{z) = Np{.z), we find from Eqs. ([7]), ([9]) that 7r(z) can also be 
written as 

oo 

= C / PM)dz'. (10) 



However, the validity of the local density approximation needs to be considered carefully. 
E.g., near the hard wall at 2; = 0, pm{z) may exhibit strong oscillations ("layering"), and also 
p{z) exhibits a nontrivial structure. Thus it is clear that the local density approximation breaks 
down near the hard wall at 2; = 0, and in fact one should use Eqs. (P). (nill) only for z > Rg, the 
gyration radius of the chains. Similarly, rapid density variations invalidating the local density 
approximation are also expected near an interface between coexisting phases (as it may occur 
for bad solvent conditions, for instance). In addition, one must require that on the scale of 
Az = Rg the change of the potential, Eq. ([6]), is negligibly small. This implies 

|A[/e,ternal(Az)/A;ijT| = XgRg/a = Rg / U < 1- (11) 



For flexible chains, Rg scales like vA^ in concentrated solutions, while for stiff chains we have 
Rg oc A^. This implies that for stiff chains considerably larger values of (and hence smaller 
values of \g) need to be chosen than for flexible ones. Of course, since one must accommodate 
in the simulation box the full density profile from a value appropriate for concentrated solutions 
near 2; = down to the dilute regime near z = H, a linear dimension H ^ C,m mandatory; 
otherwise, the distortion of the density profile due to this upper boundary at z = H leads 
to systematic errors as well. Thus, the SE method requires a careful choice of simulation 
parameters to avoid such errors and ensure the desired very good accuracy, and the large size 
H in ^-direction to some extent will reduce the advantage, that the whole equation of state can 
be estimated from a single run. 

Finally we note that the approach can be generalized to other forms of the external potential 
Uexternai{z) , different from a gravitation-like potential. E.g., if one uses 



(12) 



where k is an exponent different from one, the equation for the force balance at height z {Eq.[8]} 
gets modified as 

-{>:ja)Kp{z)z'^-\ (13) 



dp{z) 



dz 

and hence 



(^) = iCr' / '^Z''~'pm{~z)dz. (14) 



Motivation for using different potentials can be experiments made in a centrifuge. The 
centrifugal force depends linearly on the distance from the center of rotation, so that k = 2 
for this case. Although it is not clear at this point which choice for the exponent k, is optimal, 
testing that the equation of state n{pm) thus obtained does not depend on k is a nice consistency 
check. 



III. MODEL AND SIMULATION TECHNIQUE 

We use the bond fluctuation model^ on the simple cubic lattice, taking henceforth a = 1 as 
our unit of length. Each effective monomeric unit is represented by an elementary cube of the 
lattice, blocking all 8 sites at the corners of this cube from further occupation, realizing thus the 
excluded volume interaction between the monomers. The bond vectors can be taken from the 
set {(±2, 0, 0), (±2, ±1, 0), (±2, ±1, ±1), (±2, ±2, ±1), (±3, 0, 0), (±3, ±1, 0)}, including also all 
permutations between these coordinates, - altogether 108 different bond vectors occur, which 
leads to 87 different angles between successive bonds. To model the chain stiffness, an in- 
tramolecular potential depending on the angle between two successive bond vectors along the 
chain is introduced (bending energy), and also an energy term depending on the bond length 
b can be included^^, 

t^bending = + Ub = -/cOS^9(l + CCOS^9) + Eq {b - b^f . (15) 

Both the stiffness parameter / and Eq are measured in units of ksT. Note that for the model 
with / = 2.68, ^0 = 4, 6o = 0.86, c = 0.03 the isotropic to nematic transition has been studied 
by extensive Monte Carlo simulations previously^^, and in order to be able to compare the 
present results to previous work we have included this bond length energy term here again. 

Note that variable solvent quality can be modeled by including also a square well potential 
between any pair of monomers. 



fsww= r- '-^^^ (16) 




as done in Refs. ll3lJ21l . However, the present work will treat the case e = only, our analysis 
of the phase behavior when both / and e are nonzero is deferred to a later work. 

For the runs in the MVT ensemble (for the RWTI and SE methods) we define two Monte 
Carlo steps (MCS) to involve one attempt to perform a local "random hopping" move per 
each monomer unit in the system^ and one attempt of a slithering- snake move per chain. 
These are the standard moves for simulations using the bond fluctuation model^*^ in the 
canonical ensemble. The chain length used was = 20 throughout, while typical box sizes 
were 90 x 90 x 90 and 60 x 60 x 180 for the RWTI method and 80 x 80 x 250 for the SE method, 
respectively. For the SE method, a typical system contained M = 1000 chains, and about 10^ 
MCS were needed to get a reasonably well equilibrated density profile. For the RWTI method, 
the number of chains in the box was varied from very small number up to about A/" = 3000. 
Again typically 10^ MCS per run were used, taking 10^ "measurements" in each run, with 10^ 



MCS between two successive measurements (it was checked that these 10^ configurations then 
were uncorrelated) . 

For the runs in the nVT ensemble, one Monte Carlo step means one configurational bias 
move^i2ii2£i2^ plus additionally either one attempt to perform a local random hopping move 
of every effective monomer in the system or one attempt of a slithering-snake move per chain. 
For the TlfiVT method, about 50-60 different values of the chemical potential were used, and 
the box size was 90 x 90 x 90. For each parameter combination 3 ■ 10^ MCS were taken. The 
maximum value of the number of chains reached was about Af = 2700 (this corresponds to the 
polymer volume fraction of about = 0.6 in this simulation box). So the total number of MCS 
needed to record the equation of state is of the order of 2 ■ 10'' MCS (but note that 1 MCS 
needs an amount of CPU time which depends on the number of chains JV in the system). 

Since it is well knownA^^ that in the MVT ensemble the equilibration of long wavelength 
density fluctuations is very slow, suffering from "hydrodynamic slowing down"—, it is useful to 
take special precautions that the density profiles in the RWTI are well equilibrated, in order 
to avoid uncontrolled errors. Therefore for each value of the repulsive wall parameter A {cf. 
Eq. ([5])} we first equilibrated the system in the ixVT ensemble, choosing appropriately to reach 
the desired value of N . Again 3 • 10^ steps were used for this grand canonical equilibration stage. 
Then the configurational bias moves were stopped, and the MVT run with the measurement 
of the number of monomers N^aii{^) adjacent to the walls began. 

In the configurational bias moves, one needs to utilize a biased chain insertion method 
to let a polymer "grow" successively into the system. At each step all possible 108 bond 
vectors from the current effective monomer are examined, and a position for inserting the next 
monomeric unit along the chain is chosen, respecting the excluded volume condition, and using 
the Boltzmann weight calculated from the intramolecular energy, Eq. (fT5|) . The statistical 
weight of the generated polymer configuration hence is easily calculated recursively, and thus 
the bias can be accounted for in the acceptance probability for the move. 

In our simulations, we have recorded standard single-chain characteristics such as the mean- 
square end-to-end distance and the mean-square gyration radius Rg of the chains (r^ are 
positions of monomers, tcm is the position of the center of mass). 



from the runs in the ^VT ensemble. Note that these quantities depend on the polymer volume 
fraction = 8J\fNa^ /V in the system, or the average density of monomer units pm = (p/Sa^. 
The average (...) extends over all chains and all generated system configurations. In the SE 
method, where we have a density profile from a rather large density near z = to almost zero 
ai z = H, we expect that the radii Re, Rg will depend distinctly on the height zcm of the 
center of mass of a chain, and hence this information could only be obtained with substantially 
reduced statistical accuracy. In the RWTI method, one can estimate the radii as well if one 
restricts the averaging to chains with center of mass coordinate fcM sufficiently remote from 
the walls. 

IV. RESULTS FOR FULLY FLEXIBLE CHAINS 

We start by comparing results obtained from the RWTI method with results obtained from 
the thermodynamic integration in the grand canonical ensemble (Fig. [1]). The agreement be- 
tween the results of both methods actually is excellent (relative deviations are smaller than 
10~^, and the statistical error for all data points obtained by both methods was also always 
less than 1%). Also the old data by Deutsch and Dickman^, which clearly are considerably 
less accurate, agree with the present calculation to within a few percent. 




(17) 



An essential consideration for the TlfiVT method is that the values of /i for which calcu- 
lations are performed must be so closely spaced that the probability distributions P^{pm) at 
neighboring choices of /i overlap strongly. This condition has been carefully checked. Note, 
however, that the effort in CPU resources to generate the full equation of state 7r(0) in Fig. [T] 
with the TliiVT method is comparable to the effort for a single point on the equation of state 
in the RWTI method. Therefore we restrained our efforts to the accurate calculation of two 
state points with this method only. 

Fig. [2] now compares the results obtained with the SE method with those of the TlpVT 
method, and again we find perfect agreement. Note that the choice = 0.01 implies a length 

= 100 {Eq. ([7])}, while the gyration radius is Rg ~ 6.4 for = 20 and volume fractions in 
the range from 0.1 < < 0.2, which are relevant here. Thus the condition Rg -C {Eq- fHH) } 
is safely fulfilled, but nevertheless it is important to check that no visible systematic error of the 
SE method is present (to our knowledge, it is not known in which order of Rg/^m systematic 
corrections due to the gradient in density should be expected). As expected, we see some 
increase of the random statistical error with increasing 0, but the absolute magnitude of this 
error always stays clearly below 10~^, and the relative error is between 1% and 2%. 

Thus we confirm the conclusion of Addison et al.— , for a different model than studied in their 
work, that the SE method can yield a reliable estimation of the equation of state. However, care 
has to be exerted that parameters such as the height H of the simulation box and the strength 
of the external potential are chosen appropriately, and also the statistical effort needs to be 
large enough. Figs. [3] and H] contain examples of the problems that one encounters when these 
conditions are not met. E.g., when the potential is chosen too weak for the chosen size H (and 
number of chains A/"), so that the density profile (we plot here and below profiles of the polymer 
volume fraction 0) is slightly affected by the hard wall &t z = H (see Fig. [3)d), a systematic 
depression of 7r(0) in the ideal gas region (where simply = p = pm/N = (j)/8Na^ must hold) 
is found (see Fig. [3^). Conversely, when the statistics does not suffice, or equilibration time 
was too short so that the asymptotic behavior of the density profile (see Eq. [TS] below) has not 
been reached, one can also get an overestimation of the pressure in this region. Interestingly, 
even such data that are invalid in the ideal gas region still merge rather well at larger volume 
fractions, indicating the robustness of the SE method in this regime. It is also remarkable that 
the strong layering found for the potential proportional to z^^^ near the wall for A^ = 0.5 (and 
the subsequent rather rapid decrease of the density towards zero) do not create any problems, 
however. Fig. |H on the other hand, shows a case (k = 2 in the potential Eq. [T2|, A^ = 0.001) 
for which the resulting density profile (Fig. |Dd) is too steeply varying, and then it is very likely 
that the local density approximation in no longer accurate. Consequently, it is no surprise that 
the osmotic pressure comes out systematically too large (Fig. S^). However, the choice k = 2, 
X'g = 0.0001 is again in very good agreement with the result obtained for the potential Eq. ([H]) 
with Xg = 0.01 (which also agrees with the TlpVT results, as noted above). The finding that 
the three potentials, which lead to very different density profiles, nevertheless yield the same 
result for the equation of state 7r(0) is a very clear evidence that the local density approximation 
is valid, for the chosen set of parameters. 

Fig. [5] presents then a comparison of the monomer density profiles for the cases where the 
SE method yields correct results. We use here a logarithmic scale to demonstrate that for low 
densities the barometric height formula 

Pm{z) (X exp[-Uc^term.l{z) / ksT] (18) 

is fulfilled. One can see that the density profiles in Fig. [5] are compatible with Eq. (fT8l) over 
several decades in density. Only for extremely small densities the statistical scatter becomes 
significant. Thus, analyzing the data for Pm{z) in this way provides a check whether sufficient 
statistical effort has been invested. 



V. VARIABLE CHAIN STIFFNESS 



Allowing for nonzero parameters / and Eq in Eq. (1151) already has interesting effects even in 
the ideal gas limit where single chain properties dominate the behavior, and this we will discuss 
first. 

Fig. [6] shows the chemical potential per chain, fi{4>), as function of the volume fraction 
both for fully flexible and for semifiexible chains. The logarithmic scale of the abscissa is used 
to demonstrate the approach to the ideal gas law, 

/i,d(0) = ln(/) + C(/,£o), (19) 

where C{f,eo) is a constant that depends on both / and ^o- The data shown in Fig. [6] were 
obtained from simulations in the fiVT ensemble, of course. 

Fig.[7]shows the dependence of both the -bending energy, f/bending {Eq.[T5]}, and of C(/, 0) on 
the stiffness parameter /. One sees that the bending energy per chain is essentially decreasing 
linearly with / for / > 5, and then also the bending energy per monomer clearly exceeds the 
thermal energy, and hence the chains get strongly stretched (Fig.[H]). The linear variation of the 
bending energy with / in Fig. [7^ indicates that the bending energy is already fully "saturated" , 
i.e. an energetic minimum is reached. An asymptotic formula for the bending energy per chain 
for the case of fully stretched chain gives t/bending/.A/' = — 1.03(A^ — 2)/, i.e. for / = 20, iV = 20 
the bending energy per chain is equal to 370.8 (cf. Fig. [7K). 

The variation of the constant C(/, eo) with / does not have an immediately obvious in- 
terpretation, however. In order to interpret this behavior analytically, we have estimated 
^o) applying approximate single chain partition functions in terms of the independent 
trimer, quadrumer and pentamer approximations^!^. One can see (Fig. [7b) that these approx- 
imations follow the same trend as our fit to the numerical data, but clearly the convergence 
of these approximations to the numerical result from the simulation is rather slow. Obviously, 
the constant C(/, 0) reflects a delicate interplay between excluded volume and chain stiffness 
contributions to the single chain entropy that this constant measures. 

Knowing fiid{(p) we obtain fiex{4>) = ~ fJ'id{4>)y which is needed for the thermodynamic 
integration to obtain the pressure {Eqs. The result is presented in Fig. [91 One sees 

that the variation of chain stiffness at not too large volume fractions has surprisingly little effect 
on the equation of state, despite the strong change in chain extensions (Fig. [S]) and structure, 
and the decrease in the entropy of the chains (Fig. [7]). We note that the osmotic pressure for 
semifiexible chains gets slightly enhanced, in comparison to the fully flexible ones, as soon as 
one reaches a volume fraction of about ~ 0.01, where the first deviations from the ideal gas 
law TT = (j)/8Na^ = 0/160 occur (Fig. [9]d). For large 0, however, only the data for not so large 
/ lie above the curve for / = 0, while e.g. the data for / = 20 lie below those for / = if 
> 0.12. This implies that at large enough the variation of vr with / is non-monotonic: 
TT increases first, and then decreases again. Presumably this decrease reflects the onset of a 
nematic short range order — small clusters of stretched chains oriented more or less in parallel 
take less free volume than randomly oriented ones, and hence lead to a decrease of osmotic 
pressure. 

Of course, this interpretation of the pressure maximum as function of / is highly speculative, 
and a clear answer must await a careful analysis of the isotropic-nematic transition in this model. 
Figure [To] shows that our techniques are suitable to locate the transition. Here we first consider 
the case / = 2.68, Eq = 4, which was studied in our previous work on the isotropic-nematic 
transition for this model^i^. Although one expects this transition to be weakly of first order, 
and hence the 7r(0) curve should have a small two-phase coexistence region, it turns out that 
for these values of parameters the width of this two-phase coexistence region is unmeasurably 
small, on the scale of Fig. [TD] it cannot be resolved. Thus, the isotropic and nematic branches 
'^isoifJ'), T^nemilj) in Fig. [TU] meet at the transition point in the diagram almost tangentially. 



and in the 7r{(j)) curve the transition shows up only as a shghtly rounded kink. This behavior 
is compatible with previous studies of this model^^'^^. This lack of a two-phase region at the 
transition allows to carry out the integration in Eqs. from the isotropic phase over the 

transition point into the region of the nematic phase (so no reference state in the nematic phase 
is required). 

Applying the RWTI method, a very pronounced layering (extending over about 20 lattice 
units) was observed (Fig. [TTj) for the state point with the higher volume fraction (0 ^ 0.34) in 
Fig. [TOl Close to the repulsive walls, orientational ordering (or even almost crystalline packing) 
was observed. However, in the center of the system a homogeneous state at bulk density 
was clearly reached (Fig. [TTi) . and the agreement of the pressure estimates obtained (Fig. [TUil 
suggests that the observed layering (Fig. [TT]) does not invalidate the RWTI method here. 

Let us now consider the case of stiffer macromolecules, f = 8, Eq = 0. We performed grand 
canonical simulations in the cubic box L = 90, H = 90 (box 1) and elongated box L = 80, 
H = 150 (box 2). Two starting conformations have been used — the completely empty box 
and the maximally dense packed box with chains placed along one coordinate axis having bond 
lengths equal to 2. The difference between these boxes with different geometries and sizes is 
that it is possible to fill the second one (box 2) with a volume fraction equal to 1, while this is 
impossible for the first box (box 1). To characterize the orientational ordering of the bonds we 
have estimated the standard 3x3 nematic order parameter tensor 



where is the a-th component of the unit vector along the bond connecting monomers i 
and i + 1 (the largest eigenvalue of this tensor is the nematic orientational order parameter 
S). In order to distinguish between different types of nematic structures (e.g., monodomain 
vs. multidomain structures) observed in the simulatiort^i we have also calculated the largest 
eigenvalue of this tensor for each chain separately, and afterwards performed the averaging 
over all chains in the system (the single-chain orientational order parameter obtained in such a 
way is denoted S'chain)- The hysteresis for the dependencies of the density (volume fraction) 0, 
the total orientational order parameter S and of the single-chain orientational order parameter 
S'chain vs. the chemical potential for simulations in box 2 is shown in Fig. [121 

The method to locate the isotropic-nematic transition is shown in Fig. [T31 Again, the 
isotropic and nematic branches HisoifJ'), T^nemilJ^) in Fig. [T5k meet in the transition region (values 
of /i between —170 and —160) almost tangentially. Nevertheless, the intersection point of these 
two branches can be determined with a very good accuracy. The inset in Fig. [T3b shows 
the difference 7rnem(yu) — T^isoif^) in the hysteresis region. This curve crosses zero at the value 
H fa —166. The statistical error in this region was less than 0.5%, therefore we have indicated the 
1% error bars for all data points in the inset. Additionally, we have found the intersection point 
of linear fits of both branches within and in the close vicinity of the hysteresis region (their slopes 
are quite close to each other but still different). From these two procedures of data analysis we 
were able to determine the transition point as fitr = —166 it 0.5 and TTf^ = 0.026 ± 0.001. The 
value of fifr is indicated in Fig. [T21 and that of TTfr indicates the transition in Fig. [TSb where 
also the densities in coexisting isotropic and nematic phases (determined from Fig. [T2k ) are 
shown. Note, that the hysteresis in the equation of state for this value of chain stiffness is much 
broader than that in Fig. [TOl Apart from the problem, that the finite size of the box in the x- 
and y-directions may still be responsible for some systematic errors (for L = 80 the size exceeds 
a/ {RD only by about a factor 2), an accurate location of the transition and characterization of 
the discontinuities is possible. 

Now we turn to the test of the SE method for solutions of semiflexible chains (Fig. [T^ . It is 
found (see Fig. [T^ ) that for the case f = 4, Eq = the SE method is in very good agreement 
with the TlfiVT method, as in the case of flexible chains. Note however, that the regime of 
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parameters studied here does not yet encompass the nematic phase. Similar good results are 
found for the case / = 2.68, £o = 4 discussed in Fig. (TU] (to save space these data are not shown 
here) . 

However, for the case / = 8.0, Eq = some problems of SE method start to appear because 
of the formation of a nematic phase formed on the hard wall. What happens in the system is 
shown in Fig. [T5k where the profile of the orientational order parameter, S{z), is plotted together 
with the density (volume fraction) profile, 4>{z). The standard nematic order parameter tensor, 
Eq. (!20!) . was calculated for bond vectors in each of the xy-l&jeis along the z-axis separately and 
averaged over many conformations. Afterwards its maximal eigenvalue was calculated giving 
the orientational order parameter S in each layer. The points indicated by squares and circles 
in Fig. [TBh (and in the inset showing the transition regime in an enlarged scale) present the data 
for density (squares) and orientational order parameter (circles) in the isotropic (filled symbols) 
and in the nematic (open symbols) phases obtained from the grand canonical simulations at 
fi = —166 (see Fig. [12]) which exactly corresponds to the transition point. We used the values 
of the bulk densities to locate the layer z where the same value of density occurs in the system 
with the wall, and then we plotted the value of S in the bulk system at the same layer z. In 
Fig. [T5b we present the two-dimensional xz-map of the coarse-grained order parameter profile 
for this system using red, green and blue colors to represent the average local orientation of 
monomer units along x, y and z axes, respectively (the details of the calculation method can 
be found in Refj22l). From both these figures (Fig. [T^ and fTEb) one can see that the wall 
stabilizes rather thick domains of the nematically ordered phase. 

One can observe a kink in the z-dependence of the density profile exactly between the two 
coexisting bulk densities of the isotropic-nematic transition. The order parameter, however, 
starts to rise significantly before the limiting isotropic density in the bulk is reached and has 
a value around 0.5 at this density, whereas the corresponding bulk value is close to zero. 
The order parameter profile appears strongly rounded and slightly shifted with respect to the 
bulk simulations. These effects prevent an exact localization of the values of density and 
order parameter at the isotropic-nematic transition by the SE method. The rounding of the 
transition is unavoidable and mainly due to the presence of capillary wave excitations of the 
interface^, which would even increase in magnitude upon an increase of the lateral dimension 
of the simulation box— i^, in contrast to the grand canonical bulk simulations where one can 
reduce finite size effects by increasing the system size. The slight shift of the transition in the 
order parameter profile as opposed to the density profile is an indication of a precursor of a 
nematic wetting layer at the hard wall. 

It is interesting to compare the dependence of the orientational order parameter 5* on the 
density (p (see Fig. [TB]) obtained by means of SE and TlfiVT methods which can be extracted 
from the data presented in Fig. [TBk and Figs. [T^ and [T2b . respectively. Again, the TlfiVT 
data show a hysteresis while the SE data exhibit a smooth transition. It should be emphasized 
that both in the isotropic and nematic phases outside the hysteresis region the curves for both 
methods coincide with each other indicating that the SE method reproduces the properties of 
nematic phase correctly despite the vicinity of a hard wall. 

The inevitable presence of the interface leads to a smoothing of the first order transition on 
the pressure vs. density dependence (Fig. [T7k). In this Figure both the data for the SE method 
are presented using different values of H and A^, as well as the data for the Tl^VT method 
which we have already discussed above. It is clear from the principle of the SE method that it 
is impossible to observe any hysteresis on the equation of state. Fig. [T7b shows that the 7r(0) 
curves generated by the SE method start to deviate from those generated by TlfiVT method 
near = 0.28, while for < 0.28 the agreement between both methods is excellent as well 
as for (j) > 0.36 (these parts are not shown here). All curves for the SE method obtained at 
different values of H and almost coincide with each other (a variation of the lateral box linear 
dimension L has no detectable effect on SE results) except for one data set shown with small 
filled squares {H = 500, Af = 1600, A^ = 0.005). In Figs. [T7b the density profiles are presented 



and the reason for the deviation of the one data set becomes apparent: this is the only curve 
where the isotropic-nematic interface (which shows up as a characteristic kink in the density 
profile) appears very close to the hard wall (note that for / = 8.0 the end-to-end distance of 
a chain is Re ~ 40 (Fig. 8) and the chains are already rather stretched (Fig. 12b)) so that 
layering effects influence the isotropic-nematic coexistence in this case significantly. The effect 
of the layering at the wall, which can be quite strong (see Fig. [T7b) on the equation of state in 
the isotropic regime decreases with increasing distance of the isotropic-nematic interface from 
the wall, such that the curves for the SE method in (Fig. [T7k) already coincide within our 
error bars. Note also, that the ordinate of the kink in the density profiles in Fig. [T7b is the 
same for all systems at different parameters. For the bulk simulation (the TlnVT data) for 
different systems presented in Fig. [T7h we can conclude that there exists an incommensurability 
effect which influences the equation of state: large open circles show the well equilibrated data 
obtained in the cubic box L = 90, H = 90 (box 1, see also above), and the hysteresis region is 
different from the one obtained for L = 80, H = 150 (box 2) shown with stars and large open 
squares. 

Finally, we mention here a scaling of the density profiles shown in Fig. [T8l ff or details see 
the figure caption) which is also in agreement with results obtained in Ref]33l. The curves 
superimpose for the systems where the combination of parameters N'Xg has the same value. 
The curves for larger values of J^Xg look broader and smoother in these scaled variables (f)/{M\g) 
vs. z\g and show up the kink (isotropic-nematic interface) further away from the hard wall. 

It is interesting to compare our conclusion on the shift of isotropic-nematic transition with 
the results of computer simulation studies of hard rod colloidal suspensions in confinements, i.e. 
solutions confined between hard walls and/or exposed to an external gravitational potential. 
A shift of the isotropic-nematic transition to lower densities as compared to the bulk was 
found for a hard-rod fluid confined by two walls^. At the same time, a surprisingly good 
agreement between two osmotic equations of state for hard-rod fluids obtained from computer 
simulation using the SE method and from bulk simulations at many different densities has been 
reported in Ref]35l. also for densities in the nematic phase. The authors of explained this 

agreement by a very small interfacial width of the isotropic-nematic interface in comparison 
with the gravitational lengths considered in their work, a situation which is also realized in 
our simulations. However, for the model in Refj35| the density difference between isotropic and 
nematic phase is very small and any possible deviations between the bulk equation of state 
and the results from the SE method simulations in the crossover density regime are not visible 
within the resolution and statistical uncertainty of the simulations. 

Theoretically the effect of gravity on the phase behavior of hard rod solutions was studied 
Comparing density and order parameter profiles in Figs. [T^ with those calculated 



in Ref, 
in Ref. 
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49l we can see a very good agreement: the density profiles show a quite small jump 
at the transition point and a smooth but significant decrease both in nematic and isotropic 
phases, while the orientational order parameter is almost constant within each of two phases 
and experiences a quite large jump at the transition point. 



VI. CONCLUSIONS 

Monte Carlo computer simulations using the bond fluctuation model have been performed 
for solutions of semiflexible chains of the length = 20 monomer units. Three methods for 
pressure calculation in lattice Monte Carlo simulations have been investigated and compared: 
(1) the thermodynamic integration method in the grand canonical ensemble (TI/iV^T); (2) 
the repulsive wall method in the grand canonical ensemble (RWTI); (3) the sedimentation 
equilibrium method (SE) in the canonical ensemble in an external sedimentation field. All 
three methods show quite similar results for solutions of flexible chains as well as for the region 
of the isotropic phase of semiflexible chains. However, differences may occur at higher densities 



(or pressures) where for semiflexible chains the transition to the nematic phase takes place. 
Methodological problems of pressure measurement in solutions in the vicinity of isotropic- 
nematic transitions have been discussed. 

The most crucial point is that the presence of a hard repulsive wall and/or an external 
sedimentation field exerts a significant influence on the isotropic-nematic transition, both on 
the transition point (transition density) and also on the structure of the ordered phase. 

Thus we have found that the SE method is useful for obtaining the equation of state of various 
polymeric systems but its use becomes problematic in the vicinity of phase transitions. We have 
demonstrated this for the isotropic-nematic transition in solutions of semiflexible chains. The 
SE method works quite well sufficiently below and above the density of the isotropic-nematic 
transition, but it fails to predict the transition density correctly, at least for system sizes (which 
were reasonably large) and sedimentation field strengths (which were reasonably small) used in 
our simulations. 

The source of the problems that we encountered is that for a system undergoing a transition 
from the isotropic to the nematic phase the SE method implies that for densities large enough 
that the nematic phase can develop in the simulation box, one necessarily must have a transition 
zone of densities where the nematic-isotropic interface is present in the box. This nematic- 
isotropic interface is not sharp but rather extended, and hence it is not clear from data such as 
Fig. [T^ to judge where the region of the "bulk" nematic phase stops and where the region of the 
interface begins. In fact, for an equilibrium interface in the absence of an external (gravitational) 
field we would have a smooth interfacial profile between the coexisting phases as well. This 
profile can be (at least approximately) considered as the convolution of an intrinsic profile with 
capillary-wave-induced broadening, which increases with the logarithm of the lateral system 
size, proportional to InL. Note, that already in this case the problem of disentangling the 
interfacial profile from the capillary wave broadening is notoriously difficult. However, in the 
presence of a gravitational field the problem is even more subtle: while strong gravitational 
fields will lead to a "squeezing" of the intrinsic interface (similar to the squeezing of interfaces 
by very strong confinement), and eliminate the capillary wave broadening completely, weak 
gravitational fields will leave the intrinsic profile more or less intact, but eliminate the long 
wavelength part of the capillary wave spectrum. As a consequence, one expects a crossover 
from a broadening proportional to InL for not too large L to a finite width independent of 
L but controlled by the strength of the gravitational field. An explicit study of all these 
interfacial phenomena is beyond the scope of the present paper. The problem gets even more 
complicated by the fact that also in the "bulk" nematic phase the nematic order parameter 
is not at all constant, but varies with the distance z from the wall, since the order parameter 
near the transition depends strongly on density, and thus a strong variation of order parameter 
is caused by the density variation as well. In view of these problems, it is difficult from the 
SE method to estimate accurately at which density the region of the isotropic phase stops and 
at which (higher) density then the region of the nematic phase begins. An estimation of the 
nematic order parameter at the transition point is hardly possible. 

A particularly subtle difficulty occurs if parameters such as system size H in the long di- 
rection and strength of the gravitational potential are chosen such, that there is not enough 
space left for the nematic phase to develop, and the isotropic-nematic interface occurs directly 
adjacent to the wall (Fig. [T7b). Then data for the osmotic pressure in the transition region are 
obtained that are systematically too small, and suggest a transition from the isotropic to the 
nematic phase at a density that is clearly too low. This example shows that one must not rely 
on the SE method blindly when a phase transition occurs, but one needs then to check that 
the results are not changing when the strength of the gravitational field and/or the size in the 
long direction are varied. 

We restricted ourselves to the chain length = 20 monomer units because in the bond- 
fiuctuation model this chain length is sufficiently long to display Gaussian behavior in its chain 
statistics in the melt (i.e., it is a real fiexible polymer, not an oligomer) and it is not too long 



to allow for a sufficiently efficient simulation. We should emphasize that the SE method is 
valid for polymer chains of any length, provided the length scale of variation of the external 
potential in which the sedimentation equilibrium is reached is much larger than the radius of 
gyration of a polymer chain, e.g., flexible chains of the length up to 500 monomer units have 
been studied in Refj33l. However, for the semifiexible chains the effect of chain stiffening in the 
nematic state^^ will necessarily require much larger simulation boxes in comparison to the case 
of flexible chains of the same length. 

While the Tl/iVT method near the isotropic-nematic transition is hampered by hysteresis. 
Fig. [121 the fact that the transition in the (vr, fx) plane must show up as a simple intersection 
point of the curves does allow an accurate location of the transition (Fig. [T^). and to char- 
acterize the magnitudes of the jumps. In this way, the strictly horizontal part in the vr vs. 
isotherm (Fig. [T5b) can be constructed. Thus, we feel that for the accurate characterization 
of first order phase transitions the TlfiVT method is preferable, whenever applicable. With 
respect to the numerical data presented in this paper, we add the caveat that there may be 
still some systematic errors due to a too small value chosen for the lateral size L. However, the 
same size was used for the simulations by means of SE method, and hence we think that our 
discussion of the relative merits of these methods should not be affected by this problem. 

In future work, we plan to extend these studies to include also attractive interaction between 
the monomer units, to investigate the competition between nematic ordering and polymer- 
solvent phase separation, applying the methods validated in the present paper. 
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FIG. 1: Osmotic pressure tt plotted vs. polymer volume fraction (p, for the athermal fully flexible bond 
fluctuation model on the simple cubic lattice. The solid line shows the TljjbVT results, while the two 
large squares indicate results obtained with the RWTI method. Filled circles are the corresponding 
RWTI results of Deutsch and Dickman (Ref. 
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FIG. 2: Difference between osmotic pressure vr obtained with the TI/zFT and SE methods plotted vs. 
(f). Note that the difference is less than 10"'^ throughout and almost no systematic trend can be seen 
(see also discussion in Section |V] below). The SE data were obtained for a system of size 80 x 80 x 200 
containing N = 1600 chains, for a choice of \g = 0.01. 
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FIG. 3: (a) Osmotic pressure vr (on a logarithmic scale) plotted vs. volume fraction (p (on a logarithmic 
scale), for the fully athermal solution of flexible chains of length N = 20, comparing several variants of 
the SE method for a box of size 80 x 80 x 250 and M = 1200. Open squares refer to the gravitation-like 
potential, Eq. ([6]), with = 0.01, while the filled circles and stars refer to the potential in Eq. (|12p 
with K = 1/2, using = 0.1 and 0.5 correspondingly. The full line is the exactly known ideal gas 
limit, (b) Density (volume fraction) profiles of monomer units, </>(2;) = 2nm{z)/nmax, corresponding 
to the results shown in (a); solid line is for gravitation- like potential with A^ = 0.01, dashed-dotted 
line is for \'g = 0.5, k = 1/2, and dotted line is for A^ = 0.1, k = 1/2; here nm{z) is the average 
number of monomer units at the altitude z, Umax = LP' is the maximal number of monomers in 
one layer, and the coefficient 2 in the formulae for (p{z) above accounts for the fact that one fully 
occupied layer excludes for occupation all lattice sites in a neighboring layer. Note that the density 
profile for the case A^ = 0.1 is affected by the wall at z = 250 (the inset shows the same figure using a 
logarithmic scale for the density), leading to systematic errors in the equation of state. The decay of 
density profiles always occurs on scales much larger than the average value of gyration radius ^ {R'^) 
which was about 6.4 for dilute region and about 5.6 in the concentrated region close to the wall. 
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FIG. 4: (a) Same as Fig. [3^, but comparing the results for the gravitation-like potential, Eq. ([6]), with 
Xg = 0.01 (filled circles) to results for the potential in Eq. (112^ with k = 2, using = 0.0001 (open 
squares) and = 0.001 (stars), (b) Density (volume fraction) profiles of monomer units corresponding 
to the results shown in Fig. [3^; solid line is for gravitation-like potential with A^ = 0.01, dashed-dotted 
line is for A^ = 0.001, k = 2, and dotted line is for A^ = 0.0001, k = 2. 
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FIG. 5: Semi-log plot of the density profiles of monomer units for three choices of the external potential: 
Eq. Q with Xg = 0.01 (stars) with the asymptotic barometric law (solid line), Eq. ()12|) with k = 1/2, 
= 0.5 (filled circles) with the corresponding asymptotic law (dashed-dotted line) and with k = 2, 
= 0.0001 (open squares) with the asymptotic law (dashed line). For all cases the box size was 
80 X 80 X 250 and the number of chains was Af = 1200. 




,-1 



25 I . . . 1 

10"'' 10'^ 10'^ 10"^ 10° 

b) ^ 

FIG. 6: Chemical potential per chain vs. volume fraction 0, for the case of fully flexible chains, (a), 
where both / = and eo = 0, and for semiflexible ones, (b), with / = 2.68, eo = 4. Note the 
logarithmic scale of the abscissa. Straight solid line indicates the fit to Eq. (I19p . with C(0, 0) = —82.5 
and C(2.68, 4) = 38.0, respectively. 
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FIG. 7: Bending energy per chain (a) and C{f,0) (b) plotted vs. the stiffness parameter /. Thin 
soUd hne in (a) indicates the hnear fit. In part (b) C(/, 0) is plotted with stars (*) and solid line, while 

symbols indicate approximations where C{f,0) is found from a decomposition of the chain partition 
function in terms of independent trimers (open squares), tetramers (open circles) and pentamers (filled 
circles) . 
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FIG. 8: Mean squared end-to-end distance (a) and mean squared gyration radius (b) plotted vs. 
volume fraction 0, for several choices of the stiffness parameter /, / = (filled triangles), / = 1.2 
(open squares), / = 2.68 (filled squares), / = 4 (open circles), / = 8 (filled circles), / = 20 (open 
triangles), and £o = 0. 
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FIG. 9: Equation of state 7r((^) plotted vs. 0, for several choices of the stiffness parameter / and two 
choices for the bond length energy parameter eq (a). Magnified view of the equation of state in the 

low density regime (b); the ideal gas limit tt = (/)/160 is indicated by a thin solid line. The choices of 
the parameter / (and eq in part (a); part (b) refers to Eq = only) are given in the figure. Curves are 
obtained using the TlfiVT method. 
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FIG. 10: Osmotic pressure tt plotted vs. the chemical potential fi (a) and versus the volume fraction 
(f) (b), for the model with parameters / = 2.68 and £o = 4. Curves are obtained using the TliJ,VT 

method (filled triangles correspond to a dense packed starting conformation, while open squares to a 
dilute isotropic one), the two large open circles show results obtained from the RWTI method. The 
inset shows an enlarged region in the vicinity of the isotropic-nematic transition. 
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FIG. 11: Density (volume fraction) profile of a 60 x 60 x 180 system (only half of the symmetric profile 
is shown) for the parameters / = 2.68 and eo = 4 at ^ « 0.34, as used for the calculation of the 
pressure Tr{(f)) in the RWTI method. Different curves for different values of the repulsion parameter A 
from A = 0.02 to A = 1.0 are superimposed. 




FIG. 12: Hysteresis for the density vs. dependence (a), for the dependence of the single-chain 
orientational order parameter, Schaim vs. fi (b) and the total orientational order parameter, S, vs. fi 
(c). Filled triangles correspond to a dense packed starting conformation, while open squares correspond 
to a dilute isotropic starting conformation. In part (c) for the case of a dilute isotropic starting 
conformation only data points for fi below the jump to a nematic state are presented. Vertical lines 
show the estimated value of the chemical potential at the transition point, fi = —166 it 0.5 (determined 
in Fig. [13]). 
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FIG. 13: (a) The dependence of the pressure vs. chemical potential for f = 8, Eq = 0, obtained 
by TIfiVT method; the inset shows the difference TTnematic — T^isotropic in the region close to the 

isotropic-nematic transition on enlarged scales. Filled triangles correspond to a dense packed starting 
conformation, while open squares to a dilute isotropic one. (b) The equation of state for / = 8, eo = 
(solid and dotted lines). The hysteresis region as well as determined transition line are well visible. 
Two open squares indicate densities in the coexisting phases. 
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FIG. 14: (a) Equation of state vr(0) for an athermal solution of semiflexible chains, with the parameters 
of the bending energy {Eq. p^ } chosen as Eq = and / = 4 (a). Open circles are the results of the 
Tl^VT method, while the solid line is the result of the SE method, choosing the potential Eq. © and 
Xg = 0.01. The box size was equal to 80 x 80 x 250, the number of chains was equal to M = 1600. Part 
(b) shows the density profile for the SE method. The inset shows an enlarged region in the vicinity of 
the wall where the layering is well visible. 






200 



300 



400 



FIG. 15: (Color online) (a) Profiles of the orientational order parameter, S{z) (pluses +), and of the 
volume fraction of monomer units, (j){z) (crosses x), for the 80 x 80 x 1000 system, for the parameters 
/ = 8.0, Eq = 0, N = 20, M = 1600, = 0.01. Squares indicate the density values at coexistence, 
while the circles indicate the values of the nematic order parameter at coexistence. These values were 
extracted from the bulk grand canonical simulation (Fig. [T2]) . The inset shows an enlarged region with 
the isotropic-nematic interface, (b) Two dimensional xz-map of the coarse-grained order parameter 
profile for a system snapshot corresponding to (a); for details of calculation see the text and Ref 
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FIG. 16: Dependence of the nematic orientational order parameter S on the density (j): comparison of 
TlfiVT data (open and filled triangles) with SE data (solid line). The two large open circles indicate 
the isotropic-nematic transition. 
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FIG. 17: (a) Equation of state 7r((/)) for the case / = 8.0, Eq = 0, N = 
method for several choices of M, H and Xg as indicated in the legend: 



20, L = 80 using the SE 
= 0.01, H = 250, M = 1600 



(small open squares); Xg = 0.005, H = 500, M = 1600 (small filled squares); Xg = 0.005, H = 1000, 
Af = 3200 (small open circles); = 0.007, H = 500, Af = 1600 (smah filled triangles). Additionally, 
large open circles (and the solid line which is only a guide for the eye) show the result of the TlfiVT 
method for the box 90 x 90 x 90. Stars and dashed lines (as a guide for the eye) are the results of the 
TI/iVT simulations in the box 80 x 80 x 150 (from Fig. 113b). The isotropic-nematic transition line 
is visible between two large open squares indicating the coexisting densities. Regions of the density 
profiles with a kink indicating the isotropic-nematic transition are presented in (b) for A^ = 0.005 and 
different H and M as indicated in the legend. 
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FIG. 18: Rescaled density profiles (j)/{J\f\g) vs. z\g for several systems, as indicated in the legend. 
Note that curves with the same value of N\g practically superimpose. The average value of the 



gyration radius 




is about 14.8 both in the dilute and concentrated regions. 



